El termino materia obscura hace referencia a todas las especies microbianas que no han sido caracterizadas. Estas especies impiden la caracterización de ecosistemas dominados por estos microorganismos, debido a la falta de entendimiento de sus funciones.
library(igraph)
library(igraphdata)
library(networkD3)
library(RCy3)
library(tidyverse)
library(microbiomeDataSets)
library(mia)
library(phyloseq)
library(SpiecEasi)
library(ggplot2)
library(tidyr)
datos_sprockett <- SprockettTHData()
ncol(datos_sprockett)
## [1] 2086
matriz_otu <- assays(datos_sprockett)$counts %>% as.matrix()
matriz_tax <- rowData(datos_sprockett) %>% as.data.frame() %>% as.matrix()
datos_muestras <- colData(datos_sprockett) %>% as.data.frame()
OTU <- otu_table(matriz_otu, taxa_are_rows = TRUE)
TAX <- tax_table(matriz_tax)
MUESTRAS <- sample_data(datos_muestras)
physeq <- phyloseq(OTU, TAX, MUESTRAS)
physeq_especie <- tax_glom(physeq, taxrank = "Species", NArm = FALSE)
physeq_conocido <- subset_taxa(physeq_especie, !is.na(Species)) # sin materia obscura
physeq_todos <- physeq_especie # con materia obscura
filtrar_prevalencia <- function(physeq_objeto, umbral = 0.1) {
prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
conservar <- names(prevalencia[prevalencia >= umbral])
prune_taxa(conservar, physeq_objeto)
}
physeq_conocido_filtrado <- filtrar_prevalencia(physeq_conocido, 0.1)
physeq_todos_filtrado <- filtrar_prevalencia(physeq_todos, 0.1)
renombrar_especies <- function(physeq_objeto) {
tax <- tax_table(physeq_objeto)
especies <- as.character(tax[, "Species"])
especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
nombres_unicos <- make.unique(especies_limpias)
taxa_names(physeq_objeto) <- nombres_unicos
return(physeq_objeto)
}
physeq_conocido_filtrado <- renombrar_especies(physeq_conocido_filtrado)
physeq_todos_filtrado <- renombrar_especies(physeq_todos_filtrado)
red_conocido <- spiec.easi(
physeq_conocido_filtrado,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
red_todos <- spiec.easi(
physeq_todos_filtrado,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
g_conocido <- adj2igraph(getRefit(red_conocido), vertex.attr = list(name = taxa_names(physeq_conocido_filtrado)))
g_todos <- adj2igraph(getRefit(red_todos), vertex.attr = list(name = taxa_names(physeq_todos_filtrado)))
calcular_metricas <- function(grafo) {
list(
nodos = vcount(grafo),
aristas = ecount(grafo),
grado_medio = mean(degree(grafo)),
densidad = edge_density(grafo),
clustering = transitivity(grafo, type = "global"),
modularidad = modularity(cluster_louvain(grafo)),
betweenness_media = mean(betweenness(grafo))
)
}
metricas_conocido <- calcular_metricas(g_conocido)
metricas_todos <- calcular_metricas(g_todos)
tabla_comparativa <- tibble(
Métrica = names(metricas_todos),
Con_materia_obscura = unlist(metricas_todos),
Sin_materia_obscura = unlist(metricas_conocido)
)
print(tabla_comparativa)
## # A tibble: 7 × 3
## Métrica Con_materia_obscura Sin_materia_obscura
## <chr> <dbl> <dbl>
## 1 nodos 115 57
## 2 aristas 332 124
## 3 grado_medio 5.77 4.35
## 4 densidad 0.0506 0.0777
## 5 clustering 0.239 0.270
## 6 modularidad 0.516 0.533
## 7 betweenness_media 101. 45.6
plot(g_todos,
vertex.size = degree(g_todos)*2,
vertex.color = cluster_louvain(g_todos)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red con Materia Obscura")
plot(g_conocido,
vertex.size = degree(g_conocido)*2,
vertex.color = cluster_louvain(g_conocido)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red sin Materia Obscura")
eliminar_taxones_azar <- function(physeq, porcentaje = 0.3) {
taxones <- taxa_names(physeq)
eliminar <- sample(taxones, size = round(length(taxones) * porcentaje))
prune_taxa(setdiff(taxones, eliminar), physeq)
}
physeq_reducido <- eliminar_taxones_azar(physeq_conocido_filtrado)
filtrar_prevalencia <- function(physeq_objeto, umbral = 0.1) {
prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
conservar <- names(prevalencia[prevalencia >= umbral])
prune_taxa(conservar, physeq_objeto)
}
physeq_reducido_filt <- filtrar_prevalencia(physeq_reducido, 0.1)
se_reducido <- spiec.easi(
physeq_reducido_filt,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
g_reducido <- adj2igraph(getRefit(se_reducido), vertex.attr = list(name = taxa_names(physeq_reducido_filt)))
# Visualización de la red reducida
plot(g_reducido,
vertex.size = degree(g_reducido)*2,
vertex.color = cluster_louvain(g_reducido)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red con 30% de taxones eliminados")
metricas <- function(g) {
list(
nodos = vcount(g),
aristas = ecount(g),
grado_medio = mean(degree(g)),
densidad = edge_density(g),
clustering = transitivity(g, type = "global"),
modularidad = modularity(cluster_louvain(g))
)
}
m_reducido <- metricas(g_reducido)
metricas_todos <- calcular_metricas(g_todos)
metricas_conocido <- calcular_metricas(g_conocido)
metricas_bootstrap <- calcular_metricas(g_reducido)
tabla_comparativa_especie <- tibble(
Métrica = names(metricas_todos),
Con_materia_obscura = unlist(metricas_todos),
Sin_materia_obscura = unlist(metricas_conocido),
Bootstrap = unlist(metricas_bootstrap)
)
print(tabla_comparativa_especie)
## # A tibble: 7 × 4
## Métrica Con_materia_obscura Sin_materia_obscura Bootstrap
## <chr> <dbl> <dbl> <dbl>
## 1 nodos 115 57 40
## 2 aristas 332 124 54
## 3 grado_medio 5.77 4.35 2.7
## 4 densidad 0.0506 0.0777 0.0692
## 5 clustering 0.239 0.270 0.247
## 6 modularidad 0.510 0.532 0.625
## 7 betweenness_media 101. 45.6 34.6
df_metricas_especie <- tabla_comparativa_especie %>%
pivot_longer(cols = -Métrica, names_to = "Red", values_to = "Valor")
ggplot(df_metricas_especie, aes(x = Métrica, y = Valor, fill = Red)) +
geom_bar(stat = "identity", position = "dodge", color = "black") + #
scale_fill_manual(values = c("Con_materia_obscura" = "lightblue",
"Sin_materia_obscura" = "pink",
"Bootstrap" = "purple")) + # colores para cada red
theme_minimal() +
coord_flip()+
labs(title = "Comparación de métricas entre redes",
x = "Métrica",
y = "Valor de métrica") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
physeq_feces <- subset_samples(physeq, Sample_Type == "Feces")
physeq_saliva <- subset_samples(physeq, Sample_Type == "Saliva")
physeq_adultos <- subset_samples(physeq, Age_Class == "Adult")
physeq_infantes <- subset_samples(physeq, Age_Class == "Infant")
physeq_niños <- subset_samples(physeq, Age_Class == "Child")
physeq_feces <- subset_samples(physeq, Sample_Type == "Feces")
physeq_spp_f <- tax_glom(physeq_feces, taxrank = "Species", NArm = FALSE)
physeq_sindm <- subset_taxa(physeq_spp_f, !is.na(Species))
physeq_conmd <- physeq_spp_f
#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
conservar <- names(prevalencia[prevalencia >= umbral])
prune_taxa(conservar, physeq_objeto)
}
#
physeq_known_filt.feces <- prevalence_filter(physeq_sindm, 0.2)
physeq_all_filt.feces <- prevalence_filter(physeq_conmd, 0.2)
#
renombrar_especies <- function(physeq_objeto) {
tax <- tax_table(physeq_objeto)
especies <- as.character(tax[, "Species"])
# Reemplazar NA con "Desconocido"
especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
# Asegurar unicidad en los nombres
nombres_unicos <- make.unique(especies_limpias)
# Asignar los nuevos nombres como nombres de los taxones
taxa_names(physeq_objeto) <- nombres_unicos
return(physeq_objeto)
}
#
physeq_known_filt.feces <- renombrar_especies(physeq_known_filt.feces)
physeq_all_filt.feces <- renombrar_especies(physeq_all_filt.feces)
feces_known <- spiec.easi(
physeq_known_filt.feces,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
feces_all <- spiec.easi(
physeq_all_filt.feces,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
feces_known <- adj2igraph(getRefit(feces_known),
vertex.attr = list(name = taxa_names(physeq_known_filt.feces)))
feces_all <- adj2igraph(getRefit(feces_all),
vertex.attr = list(name = taxa_names(physeq_all_filt.feces)))
#
metricas <- function(g) {
list(
nodos = vcount(g),
aristas = ecount(g),
grado_medio = mean(degree(g)),
densidad = edge_density(g),
clustering = transitivity(g, type = "global"),
modularidad = modularity(cluster_louvain(g))
)
}
#
feces.m_known <- metricas(feces_known)
feces.m_all <- metricas(feces_all)
#COMPARACION:
tibble(
Métrica = names(feces.m_all),
Con_materia_obscura = unlist(feces.m_all),
Sin_materia_obscura = unlist(feces.m_known)
)
## # A tibble: 6 × 3
## Métrica Con_materia_obscura Sin_materia_obscura
## <chr> <dbl> <dbl>
## 1 nodos 51 27
## 2 aristas 92 20
## 3 grado_medio 3.61 1.48
## 4 densidad 0.0722 0.0570
## 5 clustering 0.237 0.140
## 6 modularidad 0.551 0.579
# VISUALIZACIÓN con materia obscura)
plot(feces_all,
vertex.size = degree(feces_all)*2,
vertex.color = cluster_louvain(feces_all)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red con Materia Obscura HECES")
# VISUALIZACIÓN DE g_known (sin materia obscura)
plot(feces_known,
vertex.size = degree(feces_known)*2,
vertex.color = cluster_louvain(feces_known)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red sin Materia Obscura HECES")
Red con materia oscura heces:
Red sin materia oscura heces:
physeq_saliva <- subset_samples(physeq, Sample_Type == "Saliva")
physeq_spp_s <- tax_glom(physeq_saliva, taxrank = "Species", NArm = FALSE)
physeq_saliva.sin <- subset_taxa(physeq_spp_s, !is.na(Species)) # sin materia oscura
physeq_saliva.com <- physeq_spp_s # con materia oscura
#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
conservar <- names(prevalencia[prevalencia >= umbral])
prune_taxa(conservar, physeq_objeto)
}
#
physeq_known_filt.saliva <- filtrar_prevalencia(physeq_saliva.sin, 0.2) #aplciar la funcion y sacar nuevos objetos
physeq_all_filt.saliva <- filtrar_prevalencia(physeq_saliva.com, 0.2)
#
renombrar_especies <- function(physeq_objeto) {
tax <- tax_table(physeq_objeto)
especies <- as.character(tax[, "Species"])
# Reemplazar NA con "Desconocido"
especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
# Asegurar unicidad en los nombres
nombres_unicos <- make.unique(especies_limpias)
# Asignar los nuevos nombres como nombres de los taxones
taxa_names(physeq_objeto) <- nombres_unicos
return(physeq_objeto)
}
#
physeq_known_filt.saliva <- renombrar_especies(physeq_known_filt.saliva)
physeq_all_filt.saliva <- renombrar_especies(physeq_all_filt.saliva)
# REDES SPIEC-EASI
#analsiis de redes de co-ocurrencia
saliva_known <- spiec.easi( #NO DM
physeq_known_filt.saliva,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done
saliva_all <- spiec.easi(
physeq_all_filt.saliva,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done
saliva_known <- adj2igraph(getRefit(saliva_known),
vertex.attr = list(name = taxa_names(physeq_known_filt.saliva)))
saliva_all <- adj2igraph(getRefit(saliva_all),
vertex.attr = list(name = taxa_names(physeq_all_filt.saliva)))
metricas <- function(g) {
list(
nodos = vcount(g),
aristas = ecount(g),
grado_medio = mean(degree(g)),
densidad = edge_density(g),
clustering = transitivity(g, type = "global"),
modularidad = modularity(cluster_louvain(g))
)
}
#
saliva.m_known <- metricas(saliva_known)
saliva.m_all <- metricas(saliva_all)
#COMPARACION:
tibble(
Métrica = names(saliva.m_all),
Con_materia_obscura = unlist(saliva.m_all),
Sin_materia_obscura = unlist(saliva.m_known)
)
## # A tibble: 6 × 3
## Métrica Con_materia_obscura Sin_materia_obscura
## <chr> <dbl> <dbl>
## 1 nodos 75 54
## 2 aristas 204 112
## 3 grado_medio 5.44 4.15
## 4 densidad 0.0735 0.0783
## 5 clustering 0.221 0.265
## 6 modularidad 0.484 0.520
plot(saliva_all,
vertex.size = degree(saliva_all)*2,
vertex.color = cluster_louvain(saliva_all)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red con Materia Obscura saliva")
# VISUALIZACIÓN DE g_known (sin materia obscura)
plot(saliva_known,
vertex.size = degree(saliva_known)*2,
vertex.color = cluster_louvain(saliva_known)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red sin Materia Obscura saliva")
Red con materia oscura saliva:
Red sin materia oscura saliva:
##
Redes con y sin materia obscura tomando en cuenta únicamente las
muestras de niños
physeq_niños <- subset_samples(physeq, Age_Class == "Child")
physeq_spp_niños <- tax_glom(physeq_niños, taxrank = "Species", NArm = FALSE)
physeq_niños.sin <- subset_taxa(physeq_spp_niños, !is.na(Species)) # sin materia oscura
physeq_niños.com <- physeq_spp_niños # con materia oscura
#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
conservar <- names(prevalencia[prevalencia >= umbral])
prune_taxa(conservar, physeq_objeto)
}
#
physeq_conocido_filtrado_niños <- filtrar_prevalencia(physeq_niños.sin, 0.1)
physeq_todos_filtrado_niños <- filtrar_prevalencia(physeq_niños.com, 0.1)
#
renombrar_especies <- function(physeq_objeto) {
tax <- tax_table(physeq_objeto)
especies <- as.character(tax[, "Species"])
# Reemplazar NA con "Desconocido"
especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
# Asegurar unicidad en los nombres
nombres_unicos <- make.unique(especies_limpias)
# Asignar los nuevos nombres como nombres de los taxones
taxa_names(physeq_objeto) <- nombres_unicos
return(physeq_objeto)
}
#
# Renombrar los taxa en physeq_especie (aplicable a physeq_todos y physeq_conocido)
physeq_conocido_filtrado_niños <- renombrar_especies(physeq_conocido_filtrado_niños)
physeq_todos_filtrado_niños <- renombrar_especies(physeq_todos_filtrado_niños)
# CONSTRUCCIÓN DE REDES CON SPIEC-EASI
red_conocido_niños <- spiec.easi(
physeq_conocido_filtrado_niños,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
red_todos_niños <- spiec.easi(
physeq_todos_filtrado_niños,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
g_conocido.niños <- adj2igraph(getRefit(red_conocido_niños), vertex.attr = list(name = taxa_names(physeq_conocido_filtrado_niños)))
g_todos.niños <- adj2igraph(getRefit(red_todos_niños), vertex.attr = list(name = taxa_names(physeq_todos_filtrado_niños)))
metricas <- function(g) {
list(
nodos = vcount(g),
aristas = ecount(g),
grado_medio = mean(degree(g)),
densidad = edge_density(g),
clustering = transitivity(g, type = "global"),
modularidad = modularity(cluster_louvain(g))
)
}
#
metricas_conocido_niños <- calcular_metricas(g_conocido.niños)
metricas_todos_niños <- calcular_metricas(g_todos.niños)
tabla_comparativa.niños <- tibble(
Métrica = names(metricas_todos_niños),
Con_materia_obscura = unlist(metricas_todos_niños),
Sin_materia_obscura = unlist(metricas_conocido_niños)
)
print(tabla_comparativa.niños)
## # A tibble: 7 × 3
## Métrica Con_materia_obscura Sin_materia_obscura
## <chr> <dbl> <dbl>
## 1 nodos 162 62
## 2 aristas 531 73
## 3 grado_medio 6.56 2.35
## 4 densidad 0.0407 0.0386
## 5 clustering 0.102 0.149
## 6 modularidad 0.414 0.617
## 7 betweenness_media 161. 55.1
plot(g_todos.niños,
vertex.size = degree(g_todos.niños)*2,
vertex.color = cluster_louvain(g_todos.niños)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red con Materia Obscura en niños")
plot(g_conocido.niños,
vertex.size = degree(g_conocido.niños)*2,
vertex.color = cluster_louvain(g_conocido.niños)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red sin Materia Obscura en niños")
Con materia oscura en niños (heces y saliva):
Sin materia oscura en niños (heces y saliva):
physeq_infantes <- subset_samples(physeq, Age_Class == "Infant")
physeq_spp_infantes <- tax_glom(physeq_infantes, taxrank = "Species", NArm = FALSE)
physeq_infantes.sin <- subset_taxa(physeq_spp_infantes, !is.na(Species)) # sin materia oscura
physeq_infante.com <- physeq_spp_infantes # con materia oscura
#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
conservar <- names(prevalencia[prevalencia >= umbral])
prune_taxa(conservar, physeq_objeto)
}
#
physeq_conocido_filtrado.infantes <- filtrar_prevalencia(physeq_infantes.sin, 0.1)
physeq_todos_filtrado.infantes <- filtrar_prevalencia(physeq_infante.com, 0.1)
#
renombrar_especies <- function(physeq_objeto) {
tax <- tax_table(physeq_objeto)
especies <- as.character(tax[, "Species"])
# Reemplazar NA con "Desconocido"
especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
# Asegurar unicidad en los nombres
nombres_unicos <- make.unique(especies_limpias)
# Asignar los nuevos nombres como nombres de los taxones
taxa_names(physeq_objeto) <- nombres_unicos
return(physeq_objeto)
}
#
# Renombrar los taxa en physeq_especie (aplicable a physeq_todos y physeq_conocido)
physeq_conocido_filtrado.infantes <- renombrar_especies(physeq_conocido_filtrado.infantes)
physeq_todos_filtrado.infantes <- renombrar_especies(physeq_todos_filtrado.infantes)
# CONSTRUCCIÓN DE REDES CON SPIEC-EASI
red_conocido_infantes <- spiec.easi(
physeq_conocido_filtrado.infantes,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
red_todos_infantes <- spiec.easi(
physeq_todos_filtrado.infantes,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
g_conocido.infantes <- adj2igraph(getRefit(red_conocido_infantes), vertex.attr = list(name = taxa_names(physeq_conocido_filtrado.infantes)))
g_todos.infantes <- adj2igraph(getRefit(red_todos_infantes), vertex.attr = list(name = taxa_names(physeq_todos_filtrado.infantes)))
metricas <- function(g) {
list(
nodos = vcount(g),
aristas = ecount(g),
grado_medio = mean(degree(g)),
densidad = edge_density(g),
clustering = transitivity(g, type = "global"),
modularidad = modularity(cluster_louvain(g))
)
}
#
metricas_conocido_infantes <- calcular_metricas(g_conocido.infantes)
metricas_todos_infantes <- calcular_metricas(g_todos.infantes)
tabla_comparativa.infantes <- tibble(
Métrica = names(metricas_todos_infantes),
Con_materia_obscura = unlist(metricas_todos_infantes),
Sin_materia_obscura = unlist(metricas_conocido_infantes)
)
print(tabla_comparativa.infantes)
## # A tibble: 7 × 3
## Métrica Con_materia_obscura Sin_materia_obscura
## <chr> <dbl> <dbl>
## 1 nodos 87 54
## 2 aristas 218 74
## 3 grado_medio 5.01 2.74
## 4 densidad 0.0583 0.0517
## 5 clustering 0.245 0.214
## 6 modularidad 0.504 0.518
## 7 betweenness_media 62.0 40.0
plot(g_todos.infantes,
vertex.size = degree(g_todos.infantes)*2,
vertex.color = cluster_louvain(g_todos.infantes)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red con Materia Obscura en infantes")
plot(g_conocido.infantes,
vertex.size = degree(g_conocido.infantes)*2,
vertex.color = cluster_louvain(g_conocido.infantes)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red sin Materia Obscura en infantes")
Con materia oscura en infantes (heces y saliva):
Sin materia oscura infantes (saliva y heces):
physeq_adultos <- subset_samples(physeq, Age_Class == "Adult")
physeq_spp_adultos <- tax_glom(physeq_adultos, taxrank = "Species", NArm = FALSE)
physeq_adultos.sin <- subset_taxa(physeq_spp_adultos, !is.na(Species)) # sin materia oscura
physeq_adultos.com <- physeq_spp_adultos # con materia oscura
#
prevalence_filter <- function(physeq_objeto, umbral = 0.1) {
prevalencia <- apply(otu_table(physeq_objeto), 1, function(x) mean(x > 0))
conservar <- names(prevalencia[prevalencia >= umbral])
prune_taxa(conservar, physeq_objeto)
}
#
physeq_conocido_filtrado.adultos <- filtrar_prevalencia(physeq_adultos.sin, 0.1)
physeq_todos_filtrado.adultos <- filtrar_prevalencia(physeq_adultos.com, 0.1)
#
renombrar_especies <- function(physeq_objeto) {
tax <- tax_table(physeq_objeto)
especies <- as.character(tax[, "Species"])
# Reemplazar NA con "Desconocido"
especies_limpias <- ifelse(is.na(especies), "Desconocido", especies)
# Asegurar unicidad en los nombres
nombres_unicos <- make.unique(especies_limpias)
# Asignar los nuevos nombres como nombres de los taxones
taxa_names(physeq_objeto) <- nombres_unicos
return(physeq_objeto)
}
#
# Renombrar los taxa en physeq_especie (aplicable a physeq_todos y physeq_conocido)
physeq_conocido_filtrado.adultos <- renombrar_especies(physeq_conocido_filtrado.adultos)
physeq_todos_filtrado.adultos <- renombrar_especies(physeq_todos_filtrado.adultos)
# CONSTRUCCIÓN DE REDES CON SPIEC-EASI
red_conocido_adultos <- spiec.easi(
physeq_conocido_filtrado.adultos,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done
red_todos_adultos <- spiec.easi(
physeq_todos_filtrado.adultos,
method = "mb",
lambda.min.ratio = 1e-1,
nlambda = 20,
sel.criterion = "bstars",
pulsar.params = list(thresh = 0.1)
)
## Applying data transformations...
## Selecting model with pulsar using bstars...
## Fitting final estimate with mb...
## done
g_conocido.adultos <- adj2igraph(getRefit(red_conocido_adultos),
vertex.attr = list(name = taxa_names(physeq_conocido_filtrado.adultos)))
g_todos.adultos <- adj2igraph(getRefit(red_todos_adultos),
vertex.attr = list(name = taxa_names(physeq_todos_filtrado.adultos)))
metricas <- function(g) {
list(
nodos = vcount(g),
aristas = ecount(g),
grado_medio = mean(degree(g)),
densidad = edge_density(g),
clustering = transitivity(g, type = "global"),
modularidad = modularity(cluster_louvain(g))
)
}
#
metricas_conocido_adultos <- calcular_metricas(g_conocido.adultos)
metricas_todos_adultos <- calcular_metricas(g_todos.adultos)
tabla_comparativa.adultos <- tibble(
Métrica = names(metricas_todos_adultos),
Con_materia_obscura = unlist(metricas_todos_adultos),
Sin_materia_obscura = unlist(metricas_conocido_adultos)
)
print(tabla_comparativa.adultos)
## # A tibble: 7 × 3
## Métrica Con_materia_obscura Sin_materia_obscura
## <chr> <dbl> <dbl>
## 1 nodos 177 73
## 2 aristas 737 174
## 3 grado_medio 8.33 4.77
## 4 densidad 0.0473 0.0662
## 5 clustering 0.202 0.262
## 6 modularidad 0.446 0.428
## 7 betweenness_media 164. 66.6
plot(g_todos.adultos,
vertex.size = degree(g_todos.adultos)*2,
vertex.color = cluster_louvain(g_todos.adultos)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red con Materia Obscura en adultos")
plot(g_conocido.adultos,
vertex.size = degree(g_conocido.adultos)*2,
vertex.color = cluster_louvain(g_conocido.adultos)$membership,
vertex.label.cex = 0.7,
edge.width = 1,
main = "Red sin Materia Obscura en infantes")
Con materia oscura en adultos (heces y saliva):
Sin materia oscura adultos (saliva y heces):